knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
source(here('common_fxns.R'))Create taxa-level maps of number of species impacted per cell, per year.
Identify the species included in the impacted list (using the IUCN species IDs of impact map files, from prior scripts!). Join this with the list of mapped species, to identify which taxonomic group the species is in (based on the downloaded IUCN shapefiles). Taxa groups not represented in the impacted species list will be dropped.
impact_map_dir <- file.path(dir_bd_anx, 'spp_impact_rasts')
### identify spp with impact maps
spp_impacted <- list.files(impact_map_dir, recursive = TRUE) %>%
str_extract('[0-9]+') %>%
as.integer() %>%
unique()
### Get inclusion list; drop species with no sensitivities to stressors.
### Also append the class level taxa instead of assessment group taxa
spp_sens <- get_incl_spp() %>%
filter(!is.na(stressor)) %>%
left_join(get_spp_taxon(), by = c('iucn_sid', 'assess_gp')) %>%
mutate(taxon = tolower(rank))
### look for dropped taxonomic groups...
comp_file <- here('_data', sprintf('iucn_comp_assessed_%s.csv', api_version))
spp_comp <- read_csv(comp_file)
taxa_dropped <- spp_comp %>%
filter(!assess_gp %in% spp_sens$assess_gp) %>%
.$assess_gp %>% unique()Taxonomic groups dropped from this analysis - i.e., group does not have any species that fits the “impacted” criteria:
chameleons, cycads, magnolias, fw_caridean_shrimps, fw_crayfish, conifers, surgeonfishes, crocodiles_and_alligators, fw_crabs, cacti, tarpons_and_ladyfishes, sturgeons, lobsters, amphibians
The one vulnerable tarpon species (i.e. not LC, EX, or DD) is not sensitive to any of the stressors and is not mapped, so the whole taxon is dropped from this analysis.
### directories for inputs
spp_impact_map_dir <- file.path(dir_bd_anx, 'spp_impact_rasts')
spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
### directories for outputs
taxon_impact_map_dir <- file.path(dir_bd_anx, 'taxon_impact_rasts')
### unlink(file.path(taxon_impact_map_dir, 'taxon_impact_birds_*.*'))
### unlink(file.path(taxon_impact_map_dir, '*.*'))
cell_id_rast <- raster(here('_spatial/cell_id_mol.tif'))
ocean_a_rast <- raster(here('_spatial/ocean_area_mol.tif'))These maps ignore the number of impacts occurring on each species, and count up the number of species impacted by the aggregated stressor group (land-based, ocean, climate, and fishing).
We will calculate unweighted impacts and range-weighted impacts based on protected area inclusion targets from Butchart et al 2015: 1 when range <= 1000 km2, 0.1 when range >= 250000 km2, log-linear interpolation between (no cap at 10 million km2, at this point).
For a given taxon:
get_spp_vec <- function(spp_ids, i, subgp_size) {
gp_first <- (i - 1) * subgp_size + 1
gp_last <- min(i * subgp_size, length(spp_ids))
spp_id_vec <- spp_ids[gp_first:gp_last]
}
get_spp_impact <- function(spp_id, imp_file) {
impact_map <- read_csv(imp_file, col_types = 'iii') %>%
filter(!is.na(year) & !is.na(cell_id)) %>%
filter(n_impacts > 0 & !is.na(n_impacts)) %>%
select(cell_id, year) %>%
mutate(iucn_sid = spp_id)
return(impact_map)
}
taxon_gps <- spp_sens$taxon %>% unique() %>% sort()
subgp_size <- 24 ### number of spp to process at a time
year_span <- 2003:2013
reload <- FALSEstr_cats <- c('land-based', 'ocean', 'fishing', 'climate', 'all')
### load priorities; if not normalizing, priority is based only on range, so
### easily repeatable and not dependent on the spp set
priority_df <- get_spp_priority() %>%
filter(iucn_sid %in% spp_sens$iucn_sid)
for(str_cat in str_cats) {
# str_cat <- str_cats[5]
for(taxon_gp in taxon_gps) {
### taxon_gp <- taxon_gps[6]
taxon_map_files <- file.path(taxon_impact_map_dir,
sprintf('taxon_impact_%s_%s_%s.tif',
taxon_gp, year_span, str_cat))
taxon_priority_map_files <- file.path(taxon_impact_map_dir,
sprintf('taxon_priority_impact_%s_%s_%s.tif',
taxon_gp, year_span, str_cat))
if(any(!file.exists(c(taxon_map_files, taxon_priority_map_files))) | reload) {
spp_ids_taxa <- spp_sens %>%
filter(taxon == taxon_gp) %>%
filter(str_cat == 'all' | category == str_cat) %>%
.$iucn_sid %>%
unique() %>% sort()
n_subgps <- ceiling(length(spp_ids_taxa) / subgp_size)
message('Processing ', str_cat, ' impacts for ', taxon_gp,
': ', length(spp_ids_taxa),
' species broken into ', n_subgps, ' groups.')
if(length(spp_ids_taxa) == 0) {
### uh oh, no species in this taxon affected by this stressor.
taxon_impacts_df <- data.frame()
} else {
### loop over all spp in this taxon affected by this stressor
### initialize a list to store
taxon_impacts_list <- vector('list', length = n_subgps)
for(i in 1:n_subgps) {
# i <- 1
message(' processing group ', i)
spp_id_vec <- get_spp_vec(spp_ids_taxa, i, subgp_size)
impact_maps <- parallel::mclapply(spp_id_vec,
FUN = function(spp_id) {
### spp_id <- spp_id_vec[1]
imp_file <- file.path(spp_impact_map_dir,
sprintf('spp_impact_map_%s_%s.csv', spp_id, str_cat))
imp_map <- get_spp_impact(spp_id, imp_file)
}, mc.cores = subgp_size) %>%
bind_rows()
message(' summarizing group ', i)
# system.time({
impacts_gp <- impact_maps %>%
dt_join(priority_df, by = 'iucn_sid', type = 'left') %>%
group_by(cell_id, year) %>%
summarize(priority_sum_spp = sum(priority, na.rm = TRUE),
n_spp = n()) %>%
ungroup()
# })
taxon_impacts_list[[i]] <- impacts_gp
} ### end of for loop
taxon_impacts_df <- bind_rows(taxon_impacts_list) %>%
group_by(cell_id, year) %>%
summarize(n_spp = sum(n_spp, na.rm = TRUE),
priority_sum_spp = sum(priority_sum_spp, na.rm = TRUE)) %>%
ungroup() %>%
mutate(taxon = taxon_gp)
} ### end of if statement for zero-length species vector
if(nrow(taxon_impacts_df) == 0) {
### either no species in the list, or no impacts for any of the spp
taxon_impacts_df <- data.frame(cell_id = NA,
year = NA,
n_spp = NA,
priority_sum_spp = NA,
taxon = taxon_gp)
}
### write out as count and priority rasters by year
for(i in seq_along(year_span)) {
# i <- 1
taxon_nspp_map_file <- taxon_map_files[i]
taxon_pr_map_file <- taxon_priority_map_files[i]
taxon_impact_yr_df <- taxon_impacts_df %>%
filter(year == year_span[i])
nspp_yr_rast <- map_to_rast(taxon_impact_yr_df,
cell_val = 'n_spp')
writeRaster(nspp_yr_rast, taxon_nspp_map_file, overwrite = TRUE)
pr_sum_yr_rast <- map_to_rast(taxon_impact_yr_df,
cell_val = 'priority_sum_spp')
writeRaster(pr_sum_yr_rast, taxon_pr_map_file, overwrite = TRUE)
}
} else {
message('Files exist... skipping!')
}
}
}library(animation)
make_ramp <- function(n, palette, log_out = TRUE) {
if(log_out) {
n_even <- ceiling(n)
n_max <- 10^n_even
log_steps <- c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000)
log_steps <- log_steps[log_steps <= n_max]
breaks <- seq(0, n_even, length.out = 100)
lblpos <- log10(log_steps)
labels <- 10^lblpos
colors <- hcl.colors(length(breaks), palette = palette)
} else {
n_even <- round(n + 50, -2)
breaks <- c(0, seq(.5, n_even, length.out = 100))
lblpos <- seq(0, n_even, by = 50)
labels <- lblpos
colors <- c('grey40', hcl.colors(100, palette = palette))
}
### For breaks and colors, add a different color for
### a value of exactly zero
return(list('colors' = colors,
'breaks' = breaks,
'labels' = labels,
'lblpos' = lblpos))
}
make_gifs <- function(map_stack, filename, layer_names = NULL, log_out = TRUE) {
if(is.null(layer_names)) layer_names = names(map_stack)
if(log_out) map_stack <- log10(map_stack)
n <- max(c(maxValue(map_stack), 1), na.rm = TRUE)
ramp <- make_ramp(n, palette = 'viridis', log_out)
capture.output({
saveGIF({
for(i in 1:nlayers(map_stack)){
plot(map_stack[[i]],
col = ramp$colors,
breaks = ramp$breaks,
axes = FALSE,
axis.args = list(at = ramp$lblpos, labels = ramp$labels),
main = layer_names[i])
}},
interval = 0.5, movie.name = filename,
ani.width = 700, ani.height = 420)
})
return(invisible(NULL))
}mapfiles <- list.files(taxon_impact_map_dir, full.names = TRUE)
all_taxa_maps <- data.frame(f = mapfiles) %>%
mutate(impact = str_replace_all(basename(f), 'taxon.+_[0-9]{4}_|.tif', ''),
type = ifelse(str_detect(basename(f), 'priority'), 'priority', 'sum'),
year = str_extract(basename(f), '[0-9]{4}') %>% as.integer(),
taxon = str_replace_all(basename(f), 'taxon.+impact_|_[0-9].+', ''))
str_cats <- c('fishing', 'land-based', 'ocean', 'climate', 'all')
taxa <- all_taxa_maps$taxon %>%
unique()
taxa <- taxa[c(1, 2, 3, 5, 9)]### Animate the results
for(str_cat in str_cats) {
# str_cat <- str_cats[3]
for(t in taxa) {
# t <- taxa[2]
gif_file <- here(sprintf('figs/impact_%s_%s.gif', t, str_cat))
if(!file.exists(gif_file) | reload) {
taxon_map_df <- all_taxa_maps %>%
filter(taxon == t) %>%
filter(impact == str_cat) %>%
filter(type == 'sum')
rast_files <- taxon_map_df$f
map_stack <- stack(rast_files)
make_gifs(map_stack,
filename = gif_file,
layer_names = paste(t, str_cat, 'impacts', taxon_map_df$year))
}
cat(sprintf('', gif_file))
}
}### Animate the results
for(str_cat in str_cats) {
# str_cat <- str_cats[3]
for(t in taxa) {
# t <- taxa[4]
gif_file <- here(sprintf('figs/impact_priorities_%s_%s.gif', t, str_cat))
if(!file.exists(gif_file) | reload) {
taxon_map_df <- all_taxa_maps %>%
filter(taxon == t) %>%
filter(impact == str_cat) %>%
filter(year == 2013)
rast_files <- taxon_map_df$f
priority_file <- str_detect(rast_files, 'priority')
count_rast <- raster(rast_files[!priority_file])
priority_rast <- raster(rast_files[priority_file])
### define rescaling factors
c_r <- range(values(count_rast), na.rm = TRUE)
p_r <- range(values(priority_rast), na.rm = TRUE)
if(p_r[2] - p_r[1] == 0) {
### just one value; avoid div by zero
p_rescale <- priority_rast * c_r[1] / p_r[1]
} else {
p_rescale <- (priority_rast - p_r[1]) / (p_r[2] - p_r[1])
p_rescale <- p_rescale * (c_r[2] - c_r[1]) + c_r[1]
}
map_stack <- stack(count_rast, p_rescale)
make_gifs(map_stack,
filename = gif_file,
layer_names = str_replace_all(names(map_stack), 'taxon_|.tif', ''))
}
cat(sprintf('', gif_file))
}
}